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METHOD FOR SPACE-TIME FILTERING OF NOISE IN RADIOGRAPHY 
CROSS-REFERENCE TO RELATED APPLICATIONS 

[0001] This application claims the benefit of a priority under 35 USC 1 19(a)- 

(d) to French Patent Application No. 02 13727 filed October 31, 2002, the entire 
contents of which are hereby incorporated by reference. 

BACKGROUND OF THE INVENTION 

[0002] The present invention and embodiments thereof is directed to a method 

and filter for the reduction of noise in radiography and, more particularly in 
fluoroscopic images. The invention and embodiments thereof is directed to space-time 
filtering of noise in images. The field of the invention is more particularly directed to 
the reduction of noise in images acquired in temporal sequences. 

[0003] The prior art scanner type image acquisition appeiratuses can be used to 

obtain images of an object, such as the interior of a living organism, and particularly 
the interior of a human body. Such apparatuses can be used to obtain images of 
internal organs and perform diagnoses. These images are obtained by the exposure of 
the object, such as a patient, to radiation that is received on a detector after it has gone 
through the object. The detector then produces an image that can be interpreted by a 
specialist practitioner. These images contain noise known as fluoroscopic noise. This 
fluoroscopic noise is the resultant of quantum noise, caused by the nature of the 
radiation, and electronic noise, caused by the nature of the detector. 

[0004] Fluoroscopic noise is present in the image and therefore causes 

interference with the signal, known as the useful signal or information signal, present 
in the image. This makes the interpretation of the images difficult, uncertain or even 
impossible. 

[0005] In the prior art, there are known ways of improving the signal-to-noise 

ratio, or SNR. In order to increase the ratio of information present in the image, it is 
possible to increase the intensity of the radiation. However, this makes the 
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examination more traumatic for the patient, something that is unacceptable for health 
reasons. 

10006] In the prior art, there is a known technique of filtering fluoroscopic 

noise. In this technique, a temporal mean is taken between the value of two points, or 
pixels, having the same coordinates in two images. The two images belong to one 
and the same sequence of images representing one and the same region. The two 
images have the same framing and parameters of exposure but are taken at different 
dates. If an image I is filtered at the date t, then, for each pixel of the image I, a mean 
is taken with the corresponding pixel of the image Y obtained at t-1. 

[0007] This method has several drawbacks. A first drawback lies in the low 

reduction of noise in the filtered image. This reduction is of the order of V2 at most. 
A second drawback lies in the problem of remanence, or the appearance of phantoms 
in the filtered images. Consider the acquisition of images of an artery into which it is 
known that a guide has been inserted. A guide is a cylindrical metal object that is 
introduced into the artery and is therefore visible in radiography. Given that the guide 
moves at high speed, it is possible that it will be present in the imaged zone at the date 
t-1, but not at the date t. However, since the image filtered at the date t is obtained by 
taking the average of the image acquired at the date t and the image acquired at the 
date t-1, a filtered image is obtained for the date t, and this filtered image shows a 
guide which, nevertheless, was not present at this date. Thus, an erroneous piece of 
information has been added to the image of the date t. 

[0008] A prior art solution gives a noise reduction of 15% and a remanence 

rate of 10%. This is unsatisfactory because it is a source of error or makes it 
impossible to carry out diagnosis. 

BRIEF DESCRIPTION OF THE INVENTION 

[0009] An embodiment of the invention provides for a space-time convolution 

filter to reduce noise in radiographic images. 
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[0010] An embodiment of invention is a method for filtering of noise in 

radiography comprising: 

[0011] for each pixel with coordinates (x,y) of a first image, a weighting is 

performed on the coefficients U(k,l) of a first convolution core with a dimension D, 
equivalent to a low-pass filter, as a function of a coefficient G which is a function of 
the difference computed between I(x,y) and I(x+k, y+1), where I(x,y) is the intensity 
of the pixel with coordinates (x,y) of the first image, and k and 1 are indices used to 
explore the coefficients of the convolution core, a second convolution core with 
coefficients Up(k,l) being thus obtained; 

[0012] for each pixel with coordinates (x,y) of the first image, a weighting is 

performed on the coefficients U(k,l) of the first convolution core as a function of the 
coefficient G which is a function of the difference computed between I(x,y) and 
r(x+k, y+1), where r(x,y) is the intensity of the pixel with coordinates (x,y) of a 
second image, a third convolution core with coefficients Up'(k,l) being thus obtained; 

[00131 the filtered value of I(x,y) is computed by the equation: 

F(jc, >;) = f X (r * UpikJ)J{x + ^,3; + /) + (1 - * Up\Kl)J\x ^k,y^ mi N (1) 

i = (^-X (2) 

y6[0,lj....(3) 

AT = i; z (r * + (1 - y) * up\k,i)).....i4) 

[0014] where F(x,y) is the filtered value of I(x,y). 

[0015] An embodiment of the invention is directed to a filter for carryout the 

above equation. 

[0016] An embodiment of the invention is directed to a scanner having a filter 

for carrying out the above equation or the use of such a filter in a scaimer. 
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[0017] An embodiment of the invention is directed to a computer program or a 

computer program product or an article of manufacture that is provided with a filter 
for carrying out the above equiation. 

BRIEF DESCRIPTION OF THE DRAWING 

[0018] The invention will be understood more clearly from the following 

description and from the figures that accompany it. These figures are given purely by 
way of an indication and in no way restrict the scope of the invention. Of these 
figures: 

[0019] Figure 1 is a block diagram illustrating the principles of operation of a 

filter; 

[0020] Figure 2 illustrates a frequency response function of a low-pass filter; 

[0021] Figure 3 illustrates the indexing of an image; 

[0022] Figure 4 illustrates steps of the method; and 

[0023] Figure 5 illustrates a Gaussian curve. 

DETAILED DESCRIPTION OF THE INVENTION 

[0024] Figure 1 illustrates the working of a filter in an embodiment of the 

invention. Figure 1 shows a block 101 corresponding to the implementation of the 
filter. The filter is implemented in the form of a program or specialized component. 
The fiher, or method according to an embodiment of the invention, is therefore 
implemented through instruction codes executed by a machine comprising 
microprocessor type logic circuits. These instruction codes are then recorded in a 
program memory. 

[0025] The logic components and memories then form part of a processing 

apparatus. An apparatus of this kind is either the image acquisition apparatus itself, or 
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a computer connected to the image acquisition apparatus. Here, the term "computer^* 
refers to a personal computer, workstation or the like. An image acquisition apparatus 
is, for example, a scanner. 

[0026] Figure 1 shows that the filter accepts an input of at least up to five 

parameters or inputs. In an embodiment of the invention, filter 101 may be 
considered as an apparatus implementing the method. A first input is an image 102, 
to be filtered, taken at a date t. An image of this kind may have characteristics as 
illustrated in Figure 3. Image 102 can be likened to a bitmap image. Image 102 is a 
digital image delivered by the image acquisition apparatus or obtained from signals 
delivered by the apparatus. Each pixel of the image 102 is identified by its 
coordinates (x,y) and by an intensity I(x,y). The intensity is equivalent to a grey 
level. The dynamic ranges of x,y and I depend on the detector used by the image 
acquisition apparatus. In general x varies from 0 to Xmax, and y varies from 0 to 
Ymax. Xmax and Ymax then define the resolution of the image. I(x,y) represents the 
value of the pixel having coordinates (x,y), i.e., its grey level I varies from 0 2 Imax. 

[00271 Filter 101 also accepts an input image 103 taken at the date t-1. Image 

103 belongs to the same image sequence as the image 102. In the image sequence, the 
image 103 occupies the place preceding the image 102. Images 102 and 103 have the 
same definition and represent the same region of space. In a variant embodiment of 
the invention, image 103 is the result of the filtering, of the image taken at the date t- 
1. A sequence of images is obtained during an examination by taking several 
successive shots of one and the same region, as would be done by a video camera. 
Image sequences of this kind are found when it is sought to image an organ or a 
region of an organism, during a cycle. Examples of cycles are the heart cycle, the 
respiratory cycle, or an arbitrarily fixed period of time. 

[0028] Filter 101 accepts noise statistics 104 as an input. The noise processed 

is the fluoroscopic noise for which the statistics are known as a function of the value 
of the gray level of the pixel considered. In particular, the standard deviation of the 
noise is known because it is known to be proportional to the square root of the number 
of photons reaching the detector. The number of photons is itself proportional to the 
gray level. There are known ways therefore of determining the statistics, particularly 
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the standard deviation, for fluoroscopic noise. A simple method to model a 
fluoroscopic noise present in an X-ray picture uses two successive images of one and 
the same region in such a way that the pixels of two images can be paired as a 
function of the region of space that they represent. The pairs of pixels thus obtained 
are collected together in sub-groups according to their gray level. For each sub-group, 
the standard deviation a of the values It(x, y) - It-l(x, y) is computed. A sub-group is 
distinguished by eliminating the pixels for which It(x, y) - It-l(x, y) is greater than the 
mean of the values It(x, y) - It-l(x, y) plus k times the standard deviation. The 
computations are repeated a certain number of times for each sub-group. Pairs of 
pixels (v, o) are then obtained. From these pixels, an iterative regression is performed 
to obtain a noise model according to a(v) = a.Vv + p.v + k, where v is the gray level 
and a, P and k are the coefHcients defining the noise. There are other possible 
methods such as the constitution of tables/correspondence charts used to obtain a 
statistical value for the noise. 

[0029] Filter 101 accepts firstly a parameter y and a parameter X. The 

parameter y belongs to the interval [0,1]. The lower the importance attached to the 
temporal aspect of the filter, the closer is this parameter to 1. The parameter X sets the 
strength of the filter 101. Strength X enables the tolerance of the noise to be 
modulated. The lower the value of X, the less the spatial filtering. The greater the 
value of X, the more is the spatial filtering. 

[0030] Filter 101 implements a convolution core U, or convolution core 105 

with a dimension D. A convolution curve is a square matrix comprising 
coefficients written as U(k,l) with k and 1 belonging to the interval [-(D-l)/2, ((D- 
l)/2]. When a convolution of an image I is done by the core U, it means that a new 
image F is produced for which: 

F^x,y)=f^X{Uik,l)Jix + k,y + l)) 

[0031] The core U is equivalent to a low-pass filter such as the one whose 

fi-equency response curve is shown in Figure 2. A filter such as this therefore has a 
gain of 1 for a zero firequency, and its gain decreases rapidly with the frequency. Its 
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implementation by a convolution core, in the particular case where D is equal to 5, 
obtained by taking the first coefficients of the inverse Fourier transform, gives the 
following matrix: 
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[0032] D has been chosen to be equal to 5 for reasons of computation power 

and real-time compatibility. Indeed, the filtering time must remain compatible with 
real-time use and should therefore neither put the human patient under pressure nor 
constitute a bottleneck for the image output rate. This means that it should be 
possible to process the images as fast as they are produced by the acquisition device. 
Similarly, the low-pass filter used is only one low-pass filter among the multitude of 
existing low-pass filters corresponding to an equivalent number of convolution nodes. 
The term "real-time" is understood here to mean a rate of filtering such that the 
processing rate does not exceed 100 ms per image. Ideally, this processing rate is a 
maximum of 30 ms per image. At present, this performance can be achieved for a 
reasonable cost with D equal to 5. There is no reason to rule out the possibility that in 
the future, in accordance with "Moore's law", the real-time and cost constraints can be 
met with D being equal to 7 or 9 or more. 

[0033] When the convolution core is applied to pixels located on the edge of 

the image, namely pixels such that x < (D-l)/2, or x > Xmax - (D-l)/2, or y < (D-l)/2, 
or y > Ymax - (D-l)/2, then the convolution core goes beyond the limits of the image. 
The coefficients of the convolution core corresponding to pixels located outside the 
image to be filtered are considered to correspond to pixels whose gray level value is 
arbitrarily fixed, generally at 0. 
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[0034] As an output, the filter 101 produces an image 106 corresponding to 

the filtered image 102. Image 106 is reused as the image 103 for the processing of the 
image required at the date t+1. 

[0035] Figure 4 shows a preliminary image acquisition step 401. An 

embodiment of the invention works optimally with two images belonging to one and 
the same sequence of images. However, it also works with only one image: this is the 
particular case where y equals 1. And a scanner type image acquisition apparatus 
produces image sequence. The acquisition starts at an original date that is written as 
to and is arbitrarily equal to 0. For each date tn, the image acquisition apparatus 
produces an image. All the images produced have the same definition. If not, it is 
seen to it that the image with the greatest definition has its definition reduced to that 
of the image with the smallest definition. The images can be available in digital form. 
As described here, each pixel of an image then has coordinates (x,y) and an intensity 
I(x,y) equivalent to a gray level. 

[0036] Step 401 proceeds to step 402 for processing the image t-1. In the 

present example and as embodiment, a description shall be given of the processing of 
the image acquired at the date t > tO. Step 402 therefore enables the production of the 
image 103 used by filter 101 to process the image 102 acquired at the date t. Step 402 
is identical to what will be described for the processing of the image acquired at the 
date t in the step 403. 

[0037] When the filtering is initialized, i.e., when the image acquired at the 

date tO is processed, there is no image tO-1 to play the role of the image 103. In the 
case of the first filtered image, either y is fixed at 1, or the image acquired at tO is not 
filtered. Not filtering this image amounts to considering it to be a good image and 
using it as such in the role of the image 103 for the filtering of the image acquired at 
the date tO+1. 

[0038] Step 402 proceeds to step 403 for the processing of the image I 

acquired at the date t. Step 403 comprises several sub-steps. An image is processed 
pixel by pixel. In step 404, initiating step 403, the operation therefore starts with the 
initializing of the processing: the operation takes position on the first pixel to be 
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processed. In the present example, it is considered that this is the pixel (0,0). More 
generally, description of the processing of the pixel having coordinates (x,y). Step 
404, the operation proceeds to step 405 for the computation of the values Up(k,l) for 
the pixel having coordinates (x,y). 

[0039] Step 404 is the step in which the coefficients U(k,l) of the core U are 

weighted as a function of the neighborhood of the pixel with coordinates (x,y) in the 
image I acquired at the date t. The neighborhood of a pixel with coordinates (x,y) in 
the image I is formed by the pixels of the image I closest to the pixel with 

coordinates (x,y). The pixel with coordinates (x,y) has an intensity or gray level, 
I(x,y) in the image 1. The result of this weighting is a core Up having coefficients 
Up(k,l). For the pixel with coordinates (x,y) the coefficients Up(k,l) are obtained with 
the following equation: 

Up(k,l) = U(k,l) X G(I(x+k, y+1) . I(x,y); a(I(x,y))). 

[0040] In a variant embodiment of the invention, G(e; a) is the value of the 

Gaussian curve centered on 0, with a standard deviation a and such that G(0; a) = 1. 
Figure 5 illustrates a Gaussian curve of this kind. o(I(x,y)) is the estimation of the 
standard deviation of the fluoroscopic noise for I(x,y). 

[0041] In another variant embodiment of the invention, G(8; a) is reduced to 

G(e) which is a linear function of the type G(e) = -a. e + 1, with a > 0. 

[0042] Step 405, the operation goes to step 406 in which the coefficients 

U(k,l) of the core U are weighted as a function of the neighborhood of the pixel with 
coordinates (x,y) in the image I acquired at the date t. The pixel with coordinates 
(x,y) has an intensity or gray level, r(x,y) in the image I'. The result of this weighting 
is a core Up' having coefficients Up'(k,l). For the pixel with coordinates (x,y) the 
coefficients Up'(k,l) are obtained with the following equation: 

UpXk,l) = U(k,l) X G(r(x+k, y+1) . I(x,y); X.a(I(x,y))). 

[0043] In a variant embodiment of the invention, the following equation is 

used: 

UpXk,l) = U(k,l) X G(FXx+k, y+1) - I(x,y); )i.a(I(x,y))) 
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where F*(x,y) is the intensity of the pixel having coordinates (x,y) in the image after 
filtering. 

10044] The coefficient X is used to set the strength of the filtering. In the same 

way as y, X is a parameter of the filtering in which a user can act whenever he wishes 
to do so: at the beginning of the processing of an image, during the processing of an 
image or the beginning of the processing of a sequence. There is no constraint in this 
respect. 

[0045] Step 406 the operation goes to step 407 for computing the values 

F(x,y), where F(x,y) is the intensity of the pixel having coordinates (x,y) of the image 
I after filtering. The value F(x,y) is obtained by the equation: 



[0046] In this formula, the coefGcient y is used to set the temporal dependence 

of the filter 101. y belongs to the interval [0, 1]. The closer y is to 0, the greater the 
importance attached to the image I'. At the ultimate point, if y is equal to 0, the image 
I no longer has any importance. The closer y is to 1, the lower the importance attached 
to the image W At the ultimate point, if y is equal to 1, the image Y no longer has any 
importance. If y is equal to 1 the method is equivalent to a purely spatial filter. In 
particular it can then be applied to the filtering of still images that do not belong to a 
sequence of images. 

[0047] The coefficient N is a coefficient of standardization. It is standard for 

a core and is obtained by the equation: 



[0048] In these equations L is equal to (D-l)/2, where D is the dimension of 

the core. D is therefore by nature an odd number because D-1 must be divisible by 2. 



10 



14XZ1 26398 



1 



[0049] From the step 407 the operation goes to a step 408 in which it is 

ascertained that all the pixels of the image I have truly been filtered. At the end of the 
filtering operation, there is a filtered image F for which the pixel with coordinates 
(x,y) has an intensity F(x,y), F(x,y) being a function of I(x,y) of the image I to be 
filtered. Inasmuch as the definition of the image I is known, the number of pixels to 
be processed is equal to its vertical definition by its horizontal definition. A simple 
counter can be used to define an end-of-processing criterion. It is also possible to 
process only one region of the image. In this case, the definition of the region is 
known and, therefore, the number of pixels to be processed is known. 

[0050] If there are pixels that remain to be processed, the operation goes to a 

step 409 for the iteration of x and y. If not, the operation goes to the step 410 for the 
processing of the following image. 

[0051] In step 409 the operation increases the value of x and/or of y as a 

function of the mode of scanning of the image. In practice, the image is scanned by 
lines. This means that x is increased by one unit and then, when x is equal to Xmax, y 
is increased by one unit and x is reset at 0. This process is repeated until x and y are 
equal to Xmax, and Ymax respectively. Other scanning modes are possible, for 
example scanning by columns. In general, any scanning mode enabling passage 
through all the pixels of the image/region to be filtered is acceptable for the 
implementation of the method. 

[0052] Step 409 the proceeds to step 405 for computing Up(k,l). 

[0053] At step 410 it is ascertained that the image sequence comprises an 

image t+1 to be processed/filtered. If this is the case, the operation passes from step 
410 to step 403 with the difference that the processed image is no longer the image 
acquired at the date t, but the image acquired at the date t+1, and so on and so forth 
with t+n, so long as there are images to be processed. 

[0054] Step 410 may also comprise or be followed by a step for storing the 

filtered images with a view to reusing them for printing, display or another processing 
operation. To this end, the apparatus implementing the method comprises a storage 
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memory, and/or means of comiection to a screen, and/or means of communication 
with a third party apparatus processing filtered images, and/or means for processing 
filtered images. The means for processing are, for example, instruction codes 
implemented in the same way as the method. 

[0055] The embodiments of the invention provide for the real-time 

reductionof a fluoroscopic noise. The embodiments of the invention uses a 
convolution core, equivalent to a low-pass filter, whose coefficients adapt to each 
pixel to be filtered. The adaptation is done as a function of the difference between the 
value of the pixel to be filtered and the value of the neighbors of this pixel. A 
convolution core is a square matrix that is superimposed on the image to be filtered. 
Each coefficient of the matrix then has a pixel of the image corresponding to it. The 
adaptation is also obtained as a fiinction of a known noise statistic a(v) for the value v 
of the pixel to be filtered. 

[0056] An embodiment of the invention also takes the past into consideration 

since the value of the filtered pixel is the result of the sum of two convolutions 
weighted by factors 7 and (1- y) with y in the interval [0, 1], The first convolution is 
that mentioned the previous paragraph. The second convolution is made on the image 
preceding the image to be filtered in a sequence of images. All the convolutions are 
made from the same core, weighted according to its environment. For the second 
convolution, the core is weighted as a function of the difference between the value v 
of the pixel to be filtered and the value of the pixels corresponding to the core 
superimposed on the image, and as a function ofX.o (v). 

[0057] Thus, a self-adaptive and parametrizable filtering is obtained. This is 

therefore a space-time convolution filter. The parameter y makes it possible to set the 
temporal dependence, and the parameter X makes it possible to set the strength of the 
filtering. 

[0058] One skilled in the art may make or propose various modifications to 

the structure and/or function and/or steps and/or manner and/or way and /or result and 
equivalents thereof without departing from the scope and extent of protection of the 
invention. 
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